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Abstract 

The dynamics of sheared inelastic-hard-sphere systems are studied using non-equilibrium molecular 
dynamics simulations and direct simulation Monte Carlo. In the molecular dynamics simulations Lees- 
Edwards boundary conditions are used to impose the shear. The dimensions of the simulation box are 
chosen to ensure that the systems are homogeneous and that the shear is applied uniformly. Various sys- 
tem properties are monitored, including the one-particle velocity distribution, granular temperature, stress 
tensor, collision rates, and time between collisions. The one-particle velocity distribution is found to agree 
reasonably well with an anisotropic Gaussian distribution, with only a slight overpopulation of the high ve- 
locity tails. The velocity distribution is strongly anisotropic, especially at lower densities and lower values 
of the coefficient of restitution, with the largest variance in the direction of shear. The density dependence 
of the compressibility factor of the sheared inelastic hard sphere system is quite similar to that of elastic 
hard sphere fluids. As the systems become more inelastic, the glancing collisions begin to dominate more 
direct, head-on collisions. Examination of the distribution of the time between collisions indicates that the 
collisions experienced by the particles are strongly correlated in the highly inelastic systems. A comparison 
of the simulation data is made with DSMC simulation of the Enskog equation. Results of the kinetic model 
of Montanero et al. [Montanero et al., J. Fluid Mech. 389, 391 (1999)] based on the Enskog equation are 
also included. In general, good agreement is found for high density, weakly inelastic systems. 
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I. INTRODUCTION 



In rapid granular flows [llj,|20, the mean flight time of the particles in the granular material may 
be large compared to the contact time between particles. Liter-particle interactions are modelled 
as "collisions," which play a key role in transferring momentum and other properties through the 
system. Granular materials in this flow regime can then be represented by a collection of inelastic 
hard spheres t^\4\- 

The simplicity of the inelastic hard sphere model lends itself well to theoretical analysis. In 
particular, the methods developed for the kinetic theory of equilibrium gases have been applied to 
rapidly sheared inelastic hard sphere systems. The seminal paper by Lun et al. [|5D marked the start 
of "complete" kinetic theories capable of predicting both the kinetic and coUisional properties. 
The Boltzmann equation has featured predominantly in the theory of granular gases due to its 
simpler form (e.g., see Ref J6|]); however, the most successful molecular kinetic theory to date is 
the revised Enskog theory [7], an extension to the Boltzmann equation for dense systems. Enskog 
theory assumes uncorrelated particle velocities and currently relies on a static structrual correlation 
factor from elastic fluids isj]. Approximate theories beyond the Enskog theory, such as ring theory 
iQ], have been developed and applied to granular systems, however, due to their complexity, their 
use has been limited (e.g., cooling, rare granular gases). 

Common to most kinetic theory solutions is the assumption of a steady state, spatially uniform 
distribution function. Provided scale separation exists, as is the case for elastic fluids, fluctua- 
tions from this steady state can be accounted for using the Chapman-Enskog expansion [|10|] . To 
solve the Enskog equation, approximations typically begin by taking moments of the kinetic equa- 
tion with respect to the density, velocity, and products of the velocity. These moment equations 
are used to solve for the parameters of an expansion or model. Typically, only terms up to the 
granular "temperature", or isotropic stress and rotation terms LILI. a re included as field variables. 
Anisotropic stresses can still be predicted from such a theory 1I12II : indeed, attempts have been 
made to include the full second order velocity moment fl3^ as a hydrodynamic variable to im- 
prove theoretical predictions. 

n 

Grad's method II14I1 solves Enskog theory using an expansion of the distribution function about 
a reference state. This has been applied to poly-disperse granular systems i4\ and, unlike per- 
turbative solutions, does not require assumptions on the strength of the shear. Kinetic models 
are a powerful method of generating simplified kinetic equations which retain key features of the 
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original. Montanero et al. [|l5h solved an improved Bhatnagar-Gross-Krook (BGK) kinetic model 



1711 for inelastic systems. The improved BGK model approximates the coUisional term of the 



kinetic equation using the first two velocity moments, which correspond to the coUisional stress 
and energy loss, and a general relaxation term. This leads to a simplified kinetic equation. The 
solution in the low dissipation limit is particularly attractive, as it provides estimates for the sys- 
tem properties without requiring numerical solution and compares favourably to Direct Simulation 
Monte Carlo (DSMC) resuks. 

DSMC [|l8n is a numerical simulation technique to directly solve the Boltzmann equation with- 
out requiring further approximations. This can then be used to rapidly test solutions of the kinetic 



equation. The method has a 
sheared inelastic systems [llSL 



ready been extended to the Enskog equation for homogeneously 
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m. 



While kinetic theories do offer insight into the behavior of granular materials, they are necessar- 
ily approximate. The Boltzmann and Enskog kinetic theories do not include velocity or dynamic 
structural correlations. Ring theory ^ is capable of including particle correlations; however, 
further approximations are required to make the resulting theory tractable. These correlations are 
present in moderately dense to dense systems of elastic particles, but they are enhanced by the clus- 
tering in inelastic systems [12 ll . 12211 . The failure of Boltzmann and Enskog theories at high densities 
is therefore expected, even for elastic hard sphere systems. On the other hand, non-equilibrium 
molecular dynamics (NEMD) simulations Is, 231 can, in principle, give "exact" results for driven 



inelastic-hard- sphere systems 124 . 



2511 . These simulations can be used to validate kinetic theories 



against the underlying model. Initial studies of sheared granular systems used moving bound- 
aries [13|,|26|1, such as rough walls, to introduce energy into the system. Due to the computational 
limitations, the wall separation is typically of the order of a few particle diameters, and wall ef- 
fects dominate the simulation results. For large system sizes, shear instability is observed ll27ll . 
Consequently, the results for wall driven simulations are strongly dependent on system size. 

Another manner to introduce shear in non-equilibrium molecular dynamics is the Lees-Edwards 
12811 or "sliding brick" boundary conditions. Simulations of inelastic hard-sphere systems using 



Lees-Edwards boundary conditions 123 



3211 lessen the influence of wall effects, by 



elimination of the surface of the system, but these simulations still introduce shear in an inho- 
mogenous manner, which leads to clustering instabilities [331 for larger systems. 

While there are many interesting similarities between elastic hard- sphere fluids and driven in- 
elastic hard-sphere systems, there are key differences. One is the tendency of inelastic hard spheres 
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to form clusters and patterns, while elastic hard sphere fluids tend to remain isotropic. Another ex- 
ample is the velocity distribution. The velocity of elastic hard spheres is governed by the Maxwell 
distribution, which is isotropic and Gaussian. The velocity distribution of flowing inelastic hard 
spheres is, in general anisotropic [|34ll . and can show significant deviations from the Gaussian 
distribution, especially when there is clustering. 

In this work, we examine the properties of sheared inelastic -hard-sphere systems using non- 
equilibrium event-driven molecular dynamics simulations with the SLLOD algorithm combined 
with Lees-Edwards boundary conditions. Part of the purpose of this work is to investigate, at 
a particle level, the differences between the behavior of inelastic and elastic (equilibrium) hard 
sphere systems. Another purpose of this work is to provide simulation data which can be used 
to test kinetic theory predictions for the properties of these systems. A previous study by Mon- 
tanero et al. [13511 has already compared 2D and 3D simulations of binary inelastic hard spheres 
against DSMC simulations of Enskog theory. They find good agreement over the range of mass 
ratio, size ratio and inelasticity studied; however, the clustering instability present in systems of 
large numbers of highly inelastic particles appears to limit the range of inelasticity studied. As 
mentioned previously, kinetic theories for sheared granular materials are typically developed for 
the case where the system is spatially uniform and homogeneously sheared. One of the difficulties 
with comparing the predictions of the kinetic theory with the simulation data for sheared granu- 
lar materials is the formation of clusters, which makes comparison between the two problematic. 
As a consequence, care is taken in this work to ensure that the systems remain homogeneous 
and strongly inelastic systems can be accessed. In these simulations, we investigate the collision 
statistics, such as velocity distributions, collision angles, time between collisions and mean free 
paths, of sheared inelastic hard spheres. In addition, we examine the variation of various bulk 
properties of the system, such as the viscosity, mean kinetic energy, and stress, with the packing 
fraction and coefficient of restitution of the particles. We also investigate the correlations between 
the collisions, which are neglected in most kinetic theory approaches. The remainder of this paper 
is organized as follows. In Section UIl we describe the details of the granular dynamics simula- 
tions. In Section Uni we describe the details of the DSMC simulations. In Section |IVl we present 
the results of our simulation work, including a comparison with the predictions of Enskog theory. 
Finally, a summary of the main findings is provided in Section IVl 
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II. SIMULATION DETAILS 



Non-equilibrium granular dynamics simulations were performed on systems of inelastic hard 
spheres of diameter a and mass m. The system is sheared in the y-plane in the x-direction with 



a constant strain rate of 7 using the SLLOD algorithm [|36l] . In this method, shear is applied 
through the use of the Lees-Edwards sliding brick boundary conditions |28, 3^ and the velocity 
is transformed relative to a linear velocity profile. The equations of motion are 

=Vi + Vi'^Gx (1) 



dt m 



- Vy^ae^ (2) 



where Fj is the force acting on particle i, is the position of particle i, Vj is the so called peculiar 
velocity of particle i, yi is the y-coordinate of particle i, Vy^i is the y-component of the peculiar 
velocity of particle i, is a unit vector pointing in the positive x-direction, and 7 is the strain rate. 

The peculiar velocity of a particle i is defined as the difference between its lab velocity Vj and 
the local streaming velocity (the velocity of the local streamline). For simple shear, it is given by 
the linear transformation 

Vj = Vj - yaiBa, 

The peculiar velocity is related to the dispersion of the particles from the average streamlines of 
the flow. The SLLOD equations of motion are particularly convenient as the peculiar velocity is 
naturally recovered without the need for a separate co-ordinate transformation. They allow the 



possibility of thermo staffing the system [|37l] and the study of time dependent shear flows. 

In a hard-sphere system, the spheres do not experience any force between collisions. The 
equations of motion can then be solved analytically for the trajectories of the spheres between 
collisions. The evolution of the position and peculiar velocity of particle i in the system between 
collisions is 

ri{t) =rj(to) + [vi(to) + Z/i(to)7eJ(t - to) 

=ri{to)+v,{to){t-to) 
Vi(t) =Vi(to) - Vy^i(to)j(t - to)e^ (3) 

When a particle undergoes a collision, it experiences an impulse which alters its velocity. These 
collisions are instantaneous and only occur between pairs of spheres (i.e., there are no three or 



higher body collisions). The inelasticity of the hard spheres is characterized by the coefficient of 
restitution a. This is defined through the amount of kinetic energy AE lost on collision 



Ai? = ^(l-«2)(f,,-v,,f (4) 



where Vj is the velocity of particle i immediately before collision, = Vj — Vj , and Vij = r / 1 r | 
is the unit vector pointing from the center of particle j to the center of particle i. 

Each collision preserves the total momentum of the particles involved; therefore, the change of 
velocities for a colliding pair of spheres i and j is given by 

V- =Vi - ^(1 + a) {vij ■ 

^'j + ^ (1 + «) i^ij ■ Vij ) hj (5) 

where the primes denote post collision values of the particle velocities. 

The coefficient of restitution a is, in general, a function of the relative velocity on collision. 
Viscoelastic models that incorporate this have been very successful in describing real systems 
such as steel spheres 113 811 . A common approximation in kinetic theory is to assume a constant 
coefficient of inelasticity, as this greatly simplifies the collision integrals while the basic physics 
is not significantly altered. A constant coefficient of restitution is used in this work to facilitate 
comparison against kinetic theory results. 

One concern for a system with a constant coefficient of restitution is the phenomenon of in- 
elastic collapse, where an infinite number of collisions occur between several spheres in a finite 
interval of time. Event-driven simulations will fail in the event of a single collapse event. In two 
dimensional, freely cooling, inelastic hard sphere systems undergo inelastic collapse with 
coefficient of restitution as high as 0.59. 

Inelastic collapse is rare in sheared systems ll40ll and is increasingly rare in higher dimensions; 
however, a near collapse situation can still cause a simulation to break down if the machine pre- 
cision is not sufficiently high to resolve a rapid series of successive collisions. In the simulations 
performed in this work, no partial or full collapse events were found, even for dense and highly 
inelastic systems. 

The simulation algorithm that we employ is a generalization of the standard event-driven 



molecular dynamics algorithm for hard spheres HI 
the use of the sliding brick boundary conditions 128 



(see Ref. jA'M ). The main modifications are 
and the SLLOD equations of motion. 
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Unlike the elastic hard sphere system, the inelastic hard-sphere system has no intrinsic time 
scale. The applied strain rate 7 sets the time scale of the system. Therefore, there are only two 
relevant dimensionless parameters: the density pa^ and the coefficient of restitution a. In this 
work, the density is varied from pa^ = 0.4 to 0.9, and the coefficient of restitution is varied from 
a = 0.4 to 0.9. 

Because the shear is imposed through the boundary conditions, the strain rate is only fixed at 
two points, separated by the entire height of the simulation box. In low density systems with large 
numbers of particles, clustering occurs [|27ll . This leads to a local variation of the strain rate in 
the system, and, consequently, the system will not be homogeneously sheared. At the onset of 
clustering, the size dependence of the system properties changes from the typical A^^^ scaling to 
a different behavior. To illustrate this, the mean free time of sheared inelastic-hard-spheres with 
a = 0.4 is shown in Fig. [T] These simulations were performed in a cubic simulation box where 
the system size was varied while holding the density constant. The "break" in the curves for the 
lower densities indicates the presence of cluster formation in the larger systems. The same general 
behavior manifests in all system properties and is relatively easy to detect. 

Kinetic theory studies of sheared inelastic systems typically assume that the system is homoge- 
neous and uniformly sheared, with a linear velocity profile. This makes the comparison between 
granular dynamics simulations and the kinetic theory problematic. To allow comparison with these 
theories, we ensure that the systems remain homogeneous during the course of the simulations. In 
order to avoid the clustering regime while still maintaining a large system size to provide proper 
statistics, the x-, y-, and ^-dimensions of the simulation box are set to the ratio 14.4 : 1 : 1 and 
there are a total = 7200 spheres. This ensured that the systems remained homogeneous for all 
conditions (i.e. number of particles, coefficient of restitution, and density) that were examined in 
this work. 

At the beginning of the simulations for each set of conditions, the spheres are arranged in a face 
centered cubic lattice at the appropriate density. The velocities of the spheres are initially assigned 
from a Maxwell-Boltzmann distribution. The simulations are then run for an "equilibration" period 
of 10^ collisions. Afterwards, system property data are collected over at least 10 production runs, 
each lasting 10^ collisions. The uncertainties of the data are estimated from the standard deviations 
of the results from these separate runs. In the next section, we describe the DSMC simulations 
performed. 
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FIG. 1: The system size dependence of the mean free time between coUisions tavg for simulations of a 
sheared, inelastic hard sphere system with a = 0.4 in a cubic box with sides of length Lbox- The number 
of particles is = 256 in the smallest system and N = 10, 976 in the largest. 

III. DSMC SIMULATIONS 

The DSMC method was used to numerically solve the Enskog equation. This technique has 



been described in detail previously [|15l . loll and is only covered briefly here. The peculiar velocity 
distribution function / is represented by using a collection of sample velocities or "simulated" 
particles: 



N 



/(v,t) = iV-i5^5='(v-v,(t)) 



(6) 



i=l 



where Vj(t) is the peculiar velocity of sample i at time t. At each time step At, the samples 
are evolved according to the SLLOD dynamics (see Eq. Q). The samples are then tested for 

1 (c) 

coUisional updates. At each time step, ^NFmlx pairs of samples in the collection are selected, 

(c) 

where Pmax is a parameter of the DSMC simulation. The probability that a collision between a 
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pair of samples i and j will be executed is proportional to 

Pff = 4na\p (k ■ Vi,) (k ■ V,,) At (7) 

where k is a randomly generated unit vector, = Vj — Vj — a'^kye.x is the relative lab velocity, 
6 is the Heaviside step function, and x is the radial distribution function at contact. In this work, 
the value of x is taken from the Camahan-Starling [|8|] equation of state for elastic hard spheres, 
which is given by 

X ^ ^ (8) 

(1 - u) 

where v = pna^/G is the solid fraction. 

(c) 

To optimize the simulation, the quantity Pmax is chosen to be the maximum observed value of 
P^f. This is estimated and updated during a simulation if pj^j^ exceeds Pmlx- The probability that 
a collision between samples i and j is executed is P^f* / P^l^, and, if the collision is accepted, the 
velocities are updated using Eq. (|5]) with Yij = — ak. 

For the results presented here, N = 1372 and At is selected such that ^NPmlx < 5. The dis- 
tribution functions are equilibrated for 10^ collisions, and then results are collected and averaged 
over 10 runs of 10^ collisions. 



IV. RESULTS AND DISCUSSION 



In this section, we present our simulation results for the properties of homogeneously sheared 
inelastic hard sphere systems. These results are compared against DSMC simulation of the Enskog 
equation to test the Ensko g app roximation. We also include the results from the kinetic model 



solved by Montanero et al. jlS . 



43|] . This theory is particularly interesting as it provides analytical 



results in the limit of small strain rates, along with simple expressions that approximate DSMC 
results. Without the small strain rate approximation, a more accurate numerical solution of the 
model is available lIlSll : however, the DSMC simulations already provide accurate Enskog theory 
results without further approximation. 
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A. Velocity distribution 



The kinetic energy of the system is defined through the fluctuations of the velocity of the parti- 
cles from their respective local streaming velocity 

1 ^ 

E=-Y,mvl. (9) 

k=l 

The mean kinetic energy is therefore a measure of the velocity dispersion present in the system. 
In analogy with elastic (equilibrium) hard sphere fluids, a kinetic (or "granular") temperature T is 
typically introduced through the relation 

^-NkBT={E), (10) 

where is the number of particles in the system, and fc^ is the Boltzmann constant. Although 
the physical significance of the "granular temperature" has been a subject of some controversy 



11211 . the concept has proved effective in the theoretical modelling of the properties of granular 
materials. 

The granular temperature of the sheared inelastic hard sphere system at steady state is plotted 

ecular dynamics simulations, the dotted lines 



in Fig. [21 The symbols are the results of our mo 

are the suggested expressions of Montanero et al. [|l5l] . and the solid lines are the DSMC simula- 
tion results. From the figure, it can be seen that the granular temperature of the system decreases 
with decreasing values of the coefficient of restitution. The particles in a strongly inelastic system 
rebound less from collisions; therefore, collisions in the direction of shear can quickly settle a 
particle to the velocity of the streamline. In addition, the motion of the particles off the streamline 
(in the y- and ^-directions) are more quickly dissipated by collisions with particles on neighbor- 
ing streamlines. Consequently, strongly inelastic systems have a greater tendency to follow the 
streamlines of a flow. 

At low densities, the granular temperature increases with decreasing particle densities. The 
collisions between particles transmit information regarding the mean velocity of the flow. For very 
low density systems, the collisions are relatively rare events, and between collisions a particle will 
generally travel on trajectories that deviate from the streamlines, thus contributing to the granular 
temperature. With increasing density, a particle will become increasingly "caged" by surrounding 
particles, experiencing more collisions that will keep it on a particular streamline. Therefore, one 
expects that the temperature should generally decrease with increasing particle density. However, 
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the simulation data indicate that the temperature of the system does not depend monotonically with 
the density, and a minimum is observed at a relatively high density for all the systems considered. 
The minimum becomes more pronounced as the coefficient of restitution decreases. 

We note that in dense experimental granular systems, particles mainly remain in contact with 
each other and interact by rolling or sliding past one another, rather than through collisions. In 
this regime, soft sphere models ll44ll . as opposed to hard sphere models, are more representative. 
Consequently, the applicability of the simulation results for the inelastic hard sphere system at 
high densities to experimental granular systems should be considered with care. 

In general, Enskog theory and the solution of Montanero et al. provides a fairly accurate de- 
scription of the simulation results; however, there is a large discrepancy for high values of the 
inelasticity and density. In addition, Enskog theory does not capture the presence of the minimum 
in the temperature with respect to the density. 

Equilibrium fluids obey the equipartition theorem: energy is, on average, distributed evenly 
between all degrees of freedom. In driven granular systems, however, this has been shown to not 
be the case [45|]. Figure[3^ shows the variation of {Vy)/ {vD with density, and Fig.[3h shows the 
variation of {v1)/{Vy). The symbols are the results of our molecular dynamics simulations, the 
dotted lines are the predictions of the theory of Montanero et al. [|l5|], and the solid lines are the 
DSMC simulation results. If the system obeyed the equipartition function, then both these ratios 
would be equal to one. The dispersion of the velocity parallel to the direction of shear (i.e. the 
x-direction) is consistently larger than that perpendicular to the shear, which is unsurprising as 
this is the direction in which energy is inputted to the system. The asymmetry increases with 
decreasing density and with decreasing values of the coefficient of restitution. It is interesting to 
note, however, that the fluctuations in the velocity in the y- and 2-directions are nearly equal. 

The low dissipation theory of Montanero et al. strongly under predicts the anisotropy in the 
velocity dispersion. DSMC results provide a better description but still deviate significantly from 
the simulation results at low values of the elasticity. 

Montanero et al.'s theory truncates terms within the second velocity moment of the collision 
integral and all higher terms. The full second moment could be included to improve predictions; 
however, as this is primarily a collision term it is unlikely to improve the predictions of the velocity 
anisotropy. 

The kinetic model could be expanded by relaxing to a generalized Gaussian distribution, as in 
the ellipsoidal statistical model. The extra degrees of freedom in the model would then be solved 
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FIG. 2: Variation of the mean kinetic energy per particle with density for (i) a = 0.4 (circles), (ii) a = 0.5 
(squares), (iii) a = 0.6 (diamonds), (iv) a = 0.7 (triangles-up), (v) a = 0.8 (triangles-left), and (vi) 



a = 0.9 (triangles-down). The dotted lines are the suggested expressions of Ref. 11511 and the solid lines 
are DSMC results. 

for by the inclusion of a full second velocity moment balance. This might still prove tractable and 
improve the predictions for the velocity dispersion anisotropy. 

For equilibrium systems, such as elastic -hard-sphere fluids, the velocity distribution is exactly 
given by the Maxwell-Boltzmann distribution; however, granular materials have been shown to 



deviate from this distribution [|34l 4^ |47|]. The simulation data for the distributions of the single 
particle x-, y-, and 2;-component of the peculiar velocity are shown in Fig. |4l The peculiar ve- 
locities are reduced by their mean squared values (v* = Vi/ (vf)^^^, for i = x, y, and z). The 
distributions are, in general, well described by an anisotropic Gaussian distribution. For the highly 
inelastic systems, the distributions display a slightly enhanced high velocity tail. This is most 
evident in the direction of shear (i.e. the x-direction). 
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FIG. 3: The ratio of the mean-squared- velocity (a) in the y- and x-directions and (b) in the z- and x- 
directions for sheared inelastic hard sphere sytems with (i) a = 0.4 (circles), (ii) a = 0.5 (squares), (iii) 
a = 0.6 (diamonds), (iv) a = 0.7 (triangles-up), (v) a = 0.8 (triangles-left), and (vi) a = 0.9 (triangles- 



down). The dotted lines are the suggested expressions of Ref. jlSll . and the solid lines are the DSMC 
results. 



For simulations in a cubic box at the onset of clustering in the system, the peculiar velocity 
distributions in the y- and 5;-directions can be shown to develop strong high velocity tails. In this 
case, the bulk of the particles are within a dense low strain rate zone, while the remainder reside in 
a rare, high strain rate and granular temperature region. The particles in the high strain rate region 
lead to a high velocity tail. Further studies on clustering effects are currently underway. 
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FIG. 4: The peculiar velocity distribution in the (a) x-direction, (b) y-direction, and (c) ^-direction in 
sheared inelastic hard sphere systems with (i) pa^ = 0.4 and a = 0.4 (circles), (ii) pa^ = 0.9 and a = 0.4 
(squares), (iii) pa^ = 0.4 and a = 0.9 (diamonds), and (iv) pa^ = 0.9 and a = 0.9 (triangles). The soUd 
line is a Gaussian distribution. 
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B. Stress tensor 



In this section, we examine the stress tensor. The time averaged value of the stress tensor (P) 
for a hard-sphere system is given by [|48ll 



Vt 

collisions 



N 



(11) 



Ate ^ mvk^k + (TTijiTiAvi 

k=l 

where Ate is the time interval between two consecutive collisions, Avj is the change of velocity 
of sphere i on collision, r is the time over which the stress tensor is averaged, and V is the total 
volume of the system. The first summation runs over all collisions that occur during the time r, 
and the indexes i and j refer to the spheres undergoing the collision; the index k runs over all 
particles in the system. 

The pressure p of the system, which is defined as 

is plotted in Fig. [5^. As expected, the pressure increases as the density of the system increases. 
It also increases with increasing coefficient of restitution due to the rise in granular temperature. 
The Enskog theory requires, as input, the collision rate between particles as a function of the 
density. This is typically given by the equation of state for elastic hard sphere fluids through the 
compressibility factor. The compressibility factor Z, defined as 

is plotted for the shear inelastic-hard- sphere system in Fig.|5h. The symbols represent the results 
of the simulations, and the line is the Carnahan-Starling equation of state [8] for the elastic hard 
sphere fluid. With the exception of the highest density, the compressibility factor for homoge- 
neously sheared inelastic spheres is quite similar to that for elastic hard spheres. The predictions 
of Enskog theory and Montanero et al. for the pressure (see Fig. [5^) agree fairly well with the 
simulation data. The main source of the discrepancy is due to the mis-prediction of the kinetic 
contribution to the pressure. 

The shear viscosity of a granular material is perhaps the most important design parameter in 
fast flows, quantifying the power lost per unit volume. The shear viscosity of the inelastic hard- 
sphere system was computed by two means. The first method is via the definition of the shear 
viscosity r] for simple Couette flow 

^ ^ JI^ (12) 

7 

15 



N 10 ^ 




FIG. 5: The (a) dimensionless pressure per/ (7717^) and (b) the compressibiUty factor Z for inelastic hard- 
sphere systems with (i) a = 0.4 (circles), (ii) a = 0.5 (squares), (iii) a = 0.6 (diamonds), (iv) a = 0.7 
(triangles-up), (v) a = 0.8 (triangles-left), and (vi) a = 0.9 (triangles-down). The uncertainty is smaller 



than the symbol size. The dotted lines are the suggested expressions of Ref. 11511 . and the solid lines are 
DSMC results. The solid line in (b) is the Carnahan-Starling equation of state for elastic hard sphere fluids, 
and the dotted lines are DSMC results for various a. 



An alternative method is to perform an energy balance. The work of shearing inputs energy into 
the system. Collisions between the inelastic spheres continuously dissipate kinetic energy. At 
steady-state, the average rate of energy input is equal to the average rate of energy dissipation 

m-. 



rj'j^V = -{E) (13) 
where (E) is the average rate of kinetic energy dissipation. The rate of energy dissipation is 
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directly related to the mean time between collisions tavg for a sphere by 

where is the total number of spheres in the system, and ( Ai?) is the average amount of kinetic 
energy lost per collision. 

The simulation results for the viscosity of sheared inelastic hard-sphere systems are summa- 
rized in Table m The upper entries are the values obtained from the stress tensor (see Eq. (fTTI) ). 
and the lower entries are the values obtained from the dissipation of kinetic energy (see Eq. (fT3l)). 
For all the simulation runs, the two agree within the statistical uncertainties of the simulations. 
Figure [6ti shows the dependence of the shear viscosity on the reduced density of the system, for 
various values of the coefficient of restitution. The viscosity increases with packing fraction and 
coefficient of restitution (remembering that the shear rate is equal to one). The theory of Mon- 
tanero et al. captures the full Enskog behaviour and predicts the viscosity well. Enskog theory 
deviates at low values of a and high densities where the predictions for the temperature begin to 
deviate from the simulation results (see Fig.[2l). 

In addition to the shear viscosity, we also monitor the in-plane normal stress coefficient r]- and 
the out-of-plane normal stress coefficient r/o, which are defined as Oal 



1 

2^' 



1 



The in-plane normal stress coefficient is plotted in Fig. [6b, and the out-of-plane normal stress co- 
efficient is plotted in Fig. [6]:. The simulation values deviate significantly from the predictions of 
Enskog theory; however, this is unsurprising as the velocity dispersion predictions deviate signifi- 
cantly from the simulation results (see Fig.[3l). 



C. Collision statistics 



In this section, we examine the statistics of the collisions experienced by the spheres. The mean 
time between collision tavg provides a characteristic time scale for the sheared inelastic hard sphere 
system. Figure |7] shows the variation of the mean time between collisions with the density of the 
system at various values of the coefficient of restitution. The time between collision decreases 
with a increasing density, which is expected; increasing the coefficient of restitution decreases the 
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FIG. 6: The viscosity for a homogeneously sheared inelastic-hard-sphere system with (i) a = 0.4 (circles), 
(ii) a = 0.5 (squares), (iii) a = 0.6 (diamonds), (iv) a = 0.7 (triangles-up), (v) a = 0.8 (triangles-left), 
and (vi) a = 0.9 (triangles-down). The dotted lines are the suggested expressions of Ref. jisll . and the solid 
lines are DSMC results. 

mean time between collision. The variation of tavg with the coefficient of restitution is given in the 
inset of Fig. |71 At densities roughly below pa^ = 0.6, the mean time between collision decreases 
monotonically with increasing values of the coefficient of resitution. However, at higher densities, 
there is a maximum in tavg- The Enskog theory results describe the results qualitatively well for 
low density systems but fail at high densities. 

For an elastic fluid, the velocities of different particles are, in general, uncorrected. Con- 
sequently, the velocity statistics of the individual collisions can be determined exactly. On the 
other hand, the particle velocities in a driven granular system can be strongly correlated, and their 
on-coUision statistics are not exactly known. 

The distribution of the angle 9 between the relative velocity and the relative position of two 
spheres on collision (cos^^ = rij ■ Vij / |rjj| |vij|) is given in Fig. [8l The solid line denotes an 
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FIG. 7: Mean time between collisions for sheared inelastic hard spheres with: (i) a = 0.4 (circles), (ii) 
a = 0.5 (squares), (iii) a = 0.6 (diamonds), (iv) a = 0.7 (triangles-up), (v) a = 0.8 (triangles-left), and 
(vi) a = 0.9 (triangles-down). Inset: variation of the mean time between collisions with the coefficient of 
restitution for (i) pa^ = 1.0 (circles), (ii) pa^ = 0.9 (squares), (iii) pa^ = 0.8 (diamonds), (iv) pa^ = 
0.7 (triangles-up), (v) pa^ = 0.6 (triangles-left), (vi) pa^ = 0.5 (triangles-down), and (vii) pa^ = 0.4 
(triangles-right). The solid lines are the DSMC simulation results. 



isotropic collision distribution (as is the case for elastic, elastic -hard-sphere systems). The symbols 
are the simulation data for sheared inelastic hard spheres, and the solid line is the DSMC result 
for pa^ = 0.9 and a = 0.4. For weakly inelastic systems, the distribution of the coUisional 
angle is close to that for the elastic hard-sphere system. As the inelasticity and density of the 
particles increases, however, there is a gradual increase of the frequency of "glancing" collisions 
(where cos 9 is near 0) at the expense of more "head-on" collisions (where cos 9 is close to — 1) . 
This is in agreement with the two-dimensional shearing simulation of Tan and Goldhirsch [|32n 
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and Campbell and Brennen yD. The Enskog theory does not capture this effect, as the DSMC 
simulations only display a small increased bias towards glancing collisions even in the dense, 
highly inelastic system. 

The increase in glancing collisions for strongly inelastic systems (see Fig. [8]) results primarily 
from collisions between pairs of particles orientated in the x-y plane. This occurs when the change 
of the streaming velocity over the diameter of a particle becomes significant in comparison to the 



average relative peculiar velocity [|21l1 . Particles separated in the ?/ -plane then have a significantly 
increased relative velocity which increases their probability of collision. Both the DSMC and 
granular dynamics simulation results support this; however, DSMC does not exhibit the large 
increase in collisions with a very large collision angle. 

In the inelastic hard sphere system, every collision results in a loss of kinetic energy. The 
simulation results for the distribution of the loss of kinetic energy on collision is given in Fig. [9l 
If the velocity distribution of the spheres were Gaussian (e.g., Maxwell-Boltzmann distribution), 
then the kinetic energy loss on collision would be distributed according to a Poisson distribution: 

/(AE) = -i-exp^ 



This is given by the solid line in Fig. |9l At high values of a, the distribution of the change of 
kinetic energy on collision is nearly exponential; for these systems, density does not significantly 
affect the results. 

As a decreases, there is a greater frequency of collisions that result in a very slight loss of 
kinetic energy (i.e. the initial peak in Fig. (9]). This corresponds to the increase in the glancing 
collisions in the systems. This enhancement of relatively elastic collisions is accompanied by 
an increase in collisions that result in large losses of kinetic energy (i.e. the long tail in Fig. 19]). 
These result from "head-on" collisions, which occur between particles oriented primarily in the 
x-direction where the velocity dispersion is the greatest. While these "head-on" collisions occur 
less frequently than glancing collisions in the highly inelastic systems, they are more violent. 
Increasing the density enhances these effects. 

Thus far, we have only studied the statistics of single collisions. One common assumption in 
many kinetic theories is that the individual collisions experienced by a particle are statistically 
independent. We now study the correlation between collisions by examining the time required for 
a particle to undergo a number of collisions. If the various collisions experienced by a particle can 
be considered to arrive at random times in an independent manner, then the time t required for a 
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FIG. 8: The distribution of collisional angles for (i) pa^ = 0.4 and a = 0.4 (circles), (ii) pa^ = 0.8 and 
a = 0.4 (triangles-left), (ii) pa^ = 0.9 and a = 0.4 (squares), (iii) pa^ = 0.4 and a = 0.9 (diamonds), and 
(iv) pa^ = 0.9 and a = 0.9 (triangles-up). The dashed line is for elastic hard spheres and the solid hne is 
from a DSMC simulation of pa^ = 0.9 and a = 0.4. 



particle to undergo n collisions is given by a Poisson process. The probability density function 
Pnit) that a particle experiences n collisions in a period of time t is 

(^/^avg) f t 



Pnit)= ' 7 , exp -— ) (14) 



where T{n) is the Gamma function. Deviations from this distribution are an indication of correla- 
tions between collisions. For elastic hard-sphere fluids, the Poisson process describes the collision 
time distribution fairly well, however, there are noticeable deviations, even at low densities, which 
increase with increasing density js^lsiiH]" 

The collision time distributions for homogeneously sheared inelastic hard sphere systems are 
shown in Fig. [TOl The solid lines denote the Poisson distribution, given by Eq. (fT4l) . At high 
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FIG. 9: The change of kinetic energy on collision for sheared inelastic hard sphere systems with (i) pa^ = 
0.4 and a = 0.4 (circles), (ii) pa^ = 0.9 and a = 0.4 (squares), (iii) pa^ = 0.4 and a = 0.9 (diamonds), 
and (iv) pa^ = 0.9 and a = 0.9 (triangles). The dashed Une represents the kinetic energy loss on colUsion 
if the velocity were given by the Maxwell-Boltzmann distribution, and the soUd Une is from a DSMC 
simulation at pa^ = 0.9 and a = 0.4. The inset highlights the frequency of collisions that result in low 
energy losses. 

values of the coefficient of restitution, the distributions are similar to those of elastic hard spheres 
and are fairly well described by a Poisson process. As the coefficient of restitution decreases, 
however, the simulation data deviate significantly from the Poisson process, indicating very strong 
correlations between collisions. Qualitatively, the deviations are similar to that observed for elastic 
hard sphere systems: there is an enhancement of very short and very long wait-times between 
collisions. However, these differences are much more pronounced for the inelastic hard sphere 
systems. 
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FIG. 10: Distribution of time between (a) one collision, (b) five collisions, and (c) ten collisions in sheared 

inelastic hard spheres with (i) pa^ = 0.4 and a = 0.4 (circles), (ii) pa^ = 0.9 and a = 0.4 (squares), (iii) 
pa^ = 0.4 and a = 0.9 (diamonds), and (iv) pa^ = 0.9 and a = 0.9 (triangles). The solid hne is for a 
Poisson process. 

V. CONCLUSIONS 

We have perforaied non-equilibrium molecular dynamics simulations of sheared inelastic-hard- 
systems using the SLLOD algorithm combined with Lees-Edwards boundary conditions. In these 
simulations, care was taken to ensure that the systems remain homogeneous and the shear was 
uniform across the system. As a consequence, these simulations may prove a useful reference to 
compare with the predictions of kinetic theory. 

DSMC simulations of the Enskog equation were performed to provide a solution to the ki- 
netic theory without further approximation. These compared favorably with the simulation results 
except for dense, strongly inelastic systems. The velocity anisotropy effect can be very strong 
even in homogeneous systems, and kinetic theory solutions must take this into account in their 
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approximations. 

Results were presented for the velocity statistics of individual particles in the system. The 
velocity distributions were, in general, well described by an anisotropic Gaussian. Theories based 
on the anisotropic Gaussian and the full second moment balance (e.g., see Ref. 115 3h are well 
suited to these systems. Sheared, inelastic-hard-sphere systems do not obey the equipartition 
theorem. Fluctuations of the velocity in the x-direction (the direction of shear) were greater than 
those in the y- and 5;-directions, which were both similar to each other. In addition, the granular 
temperature, which characterizes the overall fluctuation of the velocity, was observed to possess a 
minimum with respect to the density. This minimum becomes more pronounced as the coefficient 
of restitution of the spheres decreases. 

The variation of the stress in the system was also examined. The compressibility factor of the 
sheared inelastic-hard- sphere system was quite similar to that of elastic hard spheres, as estimated 
by the Camahan-Starling equation of state. The shear viscosity of the systems was computed in 
two different manners: from the average of the stress tensor and from the rate of dissipation of 
kinetic energy. The value of the viscosity from both these methods agree to within the statistical 
uncertainty of the simulations. The predictions of the Enskog equation and the kinetic theory of 
Montanero et al. [15] were in fairly good agreement with the simulation data. The in-plane and 
out-of-plane stress coefficients were also computed, but the kinetic theory predictions for these 
quantities were not as accurate. 

Finally, the collision statistics of particles in the sheared inelastic hard sphere system were 
studied. The mean time between collision was found to decrease monotonically with increasing 
density; however, at fixed density, it displays a maximum at intermediate values of the coefficient 
of restitution. Examination of the collision time distributions indicated the presence of strong cor- 
relations between collisions. Including these correlations within a kinetic theory will be important 
in developing an accurate description of high density, sheared inelastic -hard-sphere systems. 
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TABLE I: Dimensionless viscosity ria/{'m'y) of sheared inelastic hard-sphere systems at various densities 
p and coefficients of restitution a. The upper value is determined from the stress tensor (see Eq. ([TTI)). and 
the lower value is determined from the kinetic energy dissipation rate (see Eq. ([T3])). The value in brackets 
is the standard deviation over all of the runs. 



\ a 
pa^ \ 


0.4 


0.5 


0.6 


0.7 


0.8 


0.9 


0.4 


0.1054(0.0003) 


0.1294(0.0004) 


0.1625(0.0006) 


0.2134(0.0005) 


0.2970(0.0004) 


0.480(0.002) 




0.1054(0.0003) 


0.1294(0.0004) 


0.1625(0.0006) 


0.2134(0.0005) 


0.2969(0.0004) 


0.480(0.002) 


0.5 


0.1300(0.0003) 


0.1582(0.0002) 


0.1979(0.0007) 


0.2585(0.0003) 


0.361(0.001) 


0.583(0.002) 




0.1300(0.0003) 


0.1582(0.0002) 


0.1979(0.0007) 


0.2585(0.0003) 


0.361(0.001) 


0.583(0.001) 


0.6 


0.1722(0.0005) 


0.2078(0.0003) 


0.2606(0.0007) 


0.3416(0.0009) 


0.4785(0.0005) 


0.779(0.006) 




0.1722(0.0005) 


0.2078(0.0003) 


0.2606(0.0007) 


0.3417(0.0009) 


0.4785(0.0007) 


0.779(0.007) 


0.7 


0.2445(0.0004) 


0.2909(0.0008) 


0.364(0.001) 


0.4801(0.0008) 


0.677(0.002) 


1.117(0.005) 




0.2446(0.0003) 


0.2909(0.0008) 


0.364(0.001) 


0.4801(0.0009) 


0.677(0.002) 


1.117(0.005) 


0.8 


0.371(0.001) 


0.4375(0.0002) 


0.545(0.001) 


0.716(0.003) 


1.020(0.002) 


1.722(0.008) 




0.371(0.001) 


0.4375(0.0002) 


0.545(0.001) 


0.716(0.003) 


1.020(0.002) 


1.722(0.008) 


0.9 


0.643(0.007) 


0.731(0.001) 


0.885(0.002) 


1.141(0.004) 


1.64(0.01) 


2.86(0.01) 




0.643(0.008) 


0.731(0.002) 


0.885(0.002) 


1.141(0.004) 


1.64(0.01) 


2.856(0.010) 


1.0 


1.42(0.01) 


1.52(0.02) 


1.73(0.02) 


2.15(0.03) 


3.01(0.03) 


5.30(0.04) 




1.42(0.01) 


1.52(0.02) 


1.73(0.02) 


2.15(0.03) 


3.01(0.04) 


5.30(0.03) 
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